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ABSTRACT 



A quantitative evaluation of the limitations of the radiation boundary 
elements in the finite element code ATILA [Ref. 1] has been performed. Five 
three-dimensional models were employed, each representing a rigid spherical 
solid surrounded by water. Monopolar, dipolar and quadrupolar incident 
spherical waves were introduced and the corresponding scattered waves were 
computed using the ATILA code and an exact analytical solution. 

The dimensionless parameters that characterize the problem are ka, kL, 
and kR where k is the wavenumber of sound in water, a is the radius of the 
scatterer, R is the outer fluid mesh radius, and L is the thickness of the fluid 
layer. The range of values investigated were kR=1.5, 2.5, 4.0, ka=0.5, 1.0, 2.0 
and kL=0.5, 1.0. 

For axially symmetric incident fields, the maximum normalized errors 
occurred at the poles and were 9%, 12%, and 6%, respectively. Furthermore, 
the errors for monopolar and dipolar incident fields were strongly influenced by 
the location of the radiation boundary (kR), less so by the scatterer’s radius (ka); 
specifically, the error decreases with increasing kR and/or ka. The errors for 
quadrupolar incident fields do not exhibit any significant dependence on kR or 
ka. The errors for all the axially symmetric incident fields were not affected by 
variations of the element’s size (kL). For non-axially symmetric incident fields, 
the maximum deviation occurred at the equatorial points and was less than 
5.5%. 

Further investigation using a two-dimensional model is proposed in order 
to determine the range of values of ka, kL, and kR which will result in negligibly 
small errors. 
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I. INTRODUCTION 



The investigation described in this thesis is part of an ongoing research 
project in the numerical modeling of arbitrary densely packed, random volumetric 
active sonar arrays [Ref. 1]. The technique employed is an extension of the so- 
called T-matrix method, which has been applied to other scattering problems 
[Ref. 2], This method rigorously accounts for multiple scattering to all orders. To 
apply the T-matrix method to the problem of active sonar array performance 
prediction, it is necessary to compute the radiating and scattering properties of a 
single array element in a free field environment. To accomplish this for a real 
transducer, the ATILA [Ref. 3] finite element code is employed. 

ATILA is a French finite element code especially developed for the 
analysis of underwater acoustic transducers. This code was written by engineers 
at the Institut Superieur d'Electronique du Nord (ISEN), Lille, France and is in 
use by U.S. scientists within and outside the Navy working on U.S. Navy- 
sponsored research. 

Computation of the radiating and scattering properties of a sonar 
transducer using ATILA involves building a finite-element model representing the 
transducer and surrounding it by a finite-element mesh representing water, which 
is terminated by so-called “radiation boundary elements”. The radiation 
boundary elements are intended to “absorb” all incident acoustic radiation. In 
practice, they perform this function less than perfectly. This has a direct effect 
on the computation of a transducer’s radiating and scattering properties. 

The present investigation focuses on the evaluation of ATILA's radiation 
boundary elements. For this purpose, a three dimensional spherical fluid mesh, 
surrounding a spherical rigid body and terminated by radiation boundary 
elements is studied. An incident spherical wave consisting of a single 
component of a multipolar acoustic field strikes the sphere, and the resulting 
scattered wave field is calculated using ATILA. The amounts of acoustic field in 
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the intended outgoing component and in other multipolar components are then 
analyzed. The influence on the results of the fluid mesh element size, fluid mesh 
inner radius, and fluid mesh outer radius relative to the acoustic wavelength are 
investigated and quantitative guidelines developed in order to minimize the effect 
of imperfect absorption. 

The remainder of this thesis is divided into four chapters. Chapter II 
describes the theory involved in the finite element analysis, the possible types of 
analyses that can be performed, particularly the harmonic analysis of a radiating 
spherical structure excited by an incident spherical wave. 

Chapter III describes the three dimensional spherical models which were 
employed. Chapter IV presents and discusses the scattering results for different 
incident multipolar components and the influence of element size and mesh 
inner and outer radius. Chapter V concludes the thesis. Appendix A contains 
the C program used to analytically calculate the pressure scattered from the 
spherical boundary for pressure release and rigid surface boundary conditions. 
Appendix B presents the FORTRAN code used in ATI LA for the generation of 
the spherical incident pressure wave. 
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II. THEORY 



A. GENERAL PRESENTATION OF ATILA VERSION 5.03 

The finite element method is a technique that provides numerical 
solutions for boundary value problems and field calculations. [Refs. 4, 5, 6, 7, 8] 
ATILA is a finite element code developed specifically for the analysis of sonar 
transducers. 

The ATILA code is able to perform: 

1 . elastic or piezoelectric structures modeling, 

2. magnetostrictive structures modeling, 

3. periodic structures modeling, 

following a general formulation shown in Figure 1 (after [Ref. 3]), to model an 
underwater acoustic transducer. 




Figure 1. ATILA general formulation. After Ref. [3] 
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The ATILA code is based on the separation of the physical problem under 
consideration into a discrete number of elements which are in turn described by 
their nodes, in a given order. For each node there are a number of degrees of 
freedom (d.o.f.) that can be specified using certain boundary conditions. The 
elements, nodes, and the specific node-numbering order, is referred to as the 
topology of the problem. 

An ATILA job organization is carried out in several steps, as follows: 

1. Model definition. 

2. Mesh generation. 

3. Data file preparation. 

4. Running a job. 

5. Result file postprocessing. 

First of all, the type of analysis to be performed has to be specified, for 
example, harmonic analysis of a solid structure excited by an impinging wave. 
The type of elements that are required to describe the fluid and structure domain 
are then chosen. ATILA includes a library of 46 different elements to model 
composite, piezoelectric, fluid, magnetostrictive, coupling FEM-BEM, interface, 
and radiation dampers. 

The mesh generation procedure includes the assignment of coordinates 
for each node and the node-numbering order for each element. Throughout this 
procedure the whole physical problem is discretized into elements which allow 
the representation and modeling of many different geometrical shapes and lines, 
as, for example, PRIS15F, a fifteen-node isoparametric triangular base prism 
used to model homogenous fluid media, or TRIA06R, a six node isoparametric 
triangular element to prescribe a monopolar, dipolar, or multipolar radiation 
condition. 

One of the facilities of the ATILA code is the MOSAIQUE preprocessor. 
This routine enables the use of super-elements and generates all the necessary 
elements and node data for ATILA. 
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The data file includes all the neccessary input data such as the type of 
analysis, the material properties, the node coordinates, the elements, the 
boundary conditions, the loading data, and possibly the plane wave data, to carry 
out the analysis. 

Running an ATILA job involves calling up the subroutines to compute 
elementary matrices, solve equations and display the results. 

The available results file, if desired, can be postprocessed to create 
graphic displays of the structure, contours of constant value for potentials, 
pressures, and displacements. A simple flow chart of an ATILA job is shown in 
Figure 2 (after [Ref. 3]). 



1. MODEL DEFINITION 



2. MESH 



GENERATION 



3. DATA FILE PREPARATION 



<- 



4. RUNNNING A JOB 



RESTART I 
A 



5. RESULT FILE / POSTPROCESSING 



Figure 2. Flow chart of an ATILA job. After Ref.[3] 
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B. ATILA FINITE ELEMENT COMPUTATION FOR RADIATING AND 

SCATTERING PROBLEMS 

1. General Elastic or Piezoelectric Structures Modeling 

A large number of analyses can be performed. These are based upon the 
complete set of equations of elasticity in the structure, the Helmholtz equation in 
the fluid, and Poisson's equation in the elastic or piezoelectric material. 
Appropriate radiation boundary conditions are applied on the surface which 
surrounds the fluid domain. 

The unknown quantities are the nodal values of the displacement field U 
in the whole structure, the electric potential O in the piezoelectric material, and 
the pressure Pjn the fluid. The system of equations is written in matrix form: 



( [K uu ]-<a 2 [M]) [K u<[) ] -[L] 




V 




F 


[K u *r [k h ] to] 




$ 


= 


-q 


-p 2 c 2 (o 2 [L] T [0] T ([H]-co 2 [M 1 ])_ 




p 




pc 2 vj;_ 



where: 

U: vector of the nodal values of the components of the displacement field. 
O: vector of the nodal values of the electric potential. 

P: vector of the nodal values of the pressure field. 

F: vector of the nodal values of the applied forces, 
g: vector of the nodal values of the electric charges. 

\j/: vector of the nodal values of the integrated normal derivative of the 
pressure on the surface boundary S. 

[K uu ]: stiffness matrix. 

[K u<t> ]; piezoelectric matrix. 

[K h ]: dielectric matrix. 

[M]: consistent mass matrix. 

[H]: fluid (pseudo-)stiffness matrix. 
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[Mi]: consistent (pseudo-)fluid mass matrix. 

[L]: coupling matrix at the fluid structure interface (connectivity matrix). 

[0]: zero matrix, 
co: angular frequency, 
p: fluid density, 
c: fluid sound speed. 

T: means transposed. 

ATI l_A is able to perform: 

1. Static analysis of an elastic, piezoelectric, or hydroelastic system, 
where the displacement field and/or the electric potential or the pressure fields 
are required. 

2. Modal analysis of an elastic, piezoelectric, hydroelastic system, 
where the eigenfrequencies, the resonance and antiresonance frequencies, and 
the normal modes are computed. 

3. Harmonic analysis of a driven elastic or piezoelectric structure, or 
the scattering of an arbitrary incident wave by an elastic or piezoelectric 
structure. 

A scattered wave analysis of a spherical pressure wave incident on a solid 

structure is required to investigate the performance of the radiation boundary 

elements, and is presented in the following section. 

2. Harmonic Analysis of a Solid Structure Excited by an 
Impinging Wave 

In this type of analysis, the real and imaginary parts of the pressure field 
(P), the displacement field (U), the potential (O), and the electric current (icoQ) 
are computed. The pressure (P) and the flux pc 2 y, are split into incident and 

scattered parts. The normal derivative of the incident pressure on the surface 

dP 

boundary S is written ^ =[D] — where [D] is a matrix obtained by assembling the 

— 0n 

damping elements on the surface, provided by the code. 
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The incident pressure field (F^) is provided via a FORTRAN function 

"INCPRE (x,y,z,k)", included at the end of the main program by the user as 
shown in Appendix B, after [Ref 2], Utilization of a user-provided incident 
pressure allows the excitation of the structure with an incident wave of any kind, 
while the default function provided with ATILA generates a plane wave e )kx 
traveling from the positive to the negative direction along the Ox axis (e* 0 * time 
harmonic dependence is assumed). 

By assigning a specific entry in the ATILA code, we are able to compute 
the total pressure or just the scattered one, utilizing the "TOTAL" or 
"SCATTERED" attribute, respectively. 

The equations of elasticity in the structure, the Helmoltz equation in the 
fluid, and Poisson's equation in the solid material are written in matrix form : 



([K uu ]-co 2 [M]) 

[K U *] T 

-[Lf 
where: 

P es : vector of the nodal values of the elastic scattered pressure field, 

Pi : vector of the nodal values of the incident pressure field, 

[G]: frequency dependent complex linear matrix, 

[D]: matrix of the damping elements. 

Furthermore, the internal losses of the material used can be taken into 
account throughout a specified program entry "SKYLINE COMPLEX" and when 
no losses are required, with "SKYLINE REAL", respectively. 



[Ku*l 


-P-] 


[K*] 


[0] 


[Of 


[H] 

_ 2^2 , 



p C CO 



[M,] 

2 2 
pc 2 



u 

<t> 



p 2 c 2 co 2 



[G]P e 



F-PJP, 
[D] 



„ 2„2 

p CO 



[H] 



*U(- 

dn p 2 c 2 co 2 



[M| 



-)Pi 



2 2 r 1 
p 2 c 2 - 
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III. THREE DIMENSIONAL SPHERICAL MODELS 



A. INTRODUCTION 

In order to evaluate ATILA’s radiation boundary elements, we first 
developed a family of three-dimensional spherical models. These models 
simulate a spherical rigid structure surrounded by an infinite fluid environment. 

A total of five models were employed. All models are composed of an 
inner spherical boundary representing a rigid solid, several concentric spherical 
fluids layers, and an outer spherical boundary composed of radiation boundary 
elements (dampers), each representing a semi-infinite fluid region. 

B. DESCRIPTION OF MODELS 

The models developed are distinguished by their values of the inner 
radius a, the number and thickness L of each fluid layer, and the outer radius R. 
The appropriate dimensionless forms for these parameters are ka, kL, kR, where 

k = — is the acoustic wavenumber. Table 1 lists the properties of each model 
c 

used in the investigation. 

Model 1 had already been used in the calculation of the transition matrix 
for the scattering of acoustic waves from a thin elastic spherical shell [Ref. 2]. 
This model contains four fluid layers of two different thicknesses. To separately 
investigate the influence of each of the dimensionless parameters (ka, kR, kL), 
four additional models were created. 

Model 2 serves as a reference model. Each of models 3 through 5 is 
distinguished from model 2 in that the value of only one of ka, kL, kR is changed. 
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Model 1 


Model 2 


Model 3 


Model 4 


Model 5 




Original 


Basic 


ka 


kL 


kR 


Wavenumber k in meters' 1 


2 


1 


1 


1 


1 


Radius a in meters 


0.5 


0.5 


1.0 


0.5 


0.5 


Scatterers Radius ka 


1.0 


0.5 


1.0 


0.5 


0.5 


Number of Fluid Layers 


2 of 


4 of 


3 of 


2 of 


2 of 




L = 0.25 


L = 0.5 


L = 0.5 


L= 1.0 


L = 0.5 




2 of 












L = 0.5 










Element’s Size kL 


1.0 


0.5 


0.5 


1.0 


0.5 


Radiation Boundary kR 


4.0 


2.5 


2.5 


2.5 


1.5 



Table 1. Properties of each model used in the investigation. 

The solid structure is modeled by specifying a zero-flux or boundary 
condition on the scatterer surface (interface between solid and fluid domain). 
This is the default condition. 

For radiation problems, the fluid mesh outer limit must be spherical. 
Therefore, the surrounding solid structure fluid is modeled with a spherical 
surface of dimensionless radius kR = 2.5 or 1.5, for models 2 through 5. 

The fluid (water) is simulated by assuming the following properties: 

1. E=o,222*io 10 Pa (Young's modulus) 

2. p=0.1*10 4 kg/m 3 (Density) 

3. v=0.0 (Poisson’s Ratio) 

The whole solid structure and fluid mesh was constructed using fifteen- 
node, three-dimensional isoparametric right triangular prismatic elements 
(PRIS15) from the ATILA finite element library [Ref.3], 

ATILA offers monopolar, dipolar, and multipolar radiation boundary 
elements. Multipolar damping elements were selected to terminate the fluid 
mesh outer radius and the appropriate condition was set in the data input file. 
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The multipolar damping element used is a six-node isoparametric triangular 
element (TRIA06R) and it has to be spherical and attached to the outer surface 
of the fluid mesh. 

The topology of the prismatic (PRIS15) and of the triangular (TRIA06R) 
elements are shown in Figures 3 and 4, respectively. The numbers represent 
the nodes and the order of numbering required for each element. 

C. MESH GENERATION 

The major limitation of the three-dimensional finite element model is the 
number of degrees of freedom (D.O.F.) available (displacement along the 
coordinate axes, rotation along the coordinate axes, pressure, electric or 
magnetic potential, and magnetic excitation currents). For this reason the 
models developed were limited by the number of nodes and elements that were 
allowed. 

Specifically, the number of nodes employed in the models ranged from 
1546 to 2346 for the structure and fluid mesh, and the total number of elements 
ranged from 216 to 360. Of these, 144 to 288 are fluid elements and 72 are 
damping elements. The fluid mesh is arranged in successive layers of various 
thicknesses. The original fluid mesh was divided into four layers, two of them 
with 0.25m thickness and the other two fluid layers with 0.5m thickness. 

Figure 5 shows a Mercator projection of the nodes and elements on the 
inner surface of the fluid. The node numbers and selected element numbers 
corresponding to this layer are given. 

For the mesh generation procedure all the nodes were specified using 
spherical coordinates with the origin at the center of the spherical structure. The 
mesh spacing was always less than a quarter of the acoustic wavelength used. 

The aspect ratio of each element according to the reference manual 
should be less than or equal to 3; in our models it is 1 to 4. 
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Figure 3. PRIS15F Element Topology. After Ref.[3] 




Figure 4. TRIA06R Element Topology. After Ref. [3] 
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Figure 5. Inner Surface Fluid Layer Model, Mercator Projection. 
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The internal angles of our elements were modeled between 30 and 90 
degrees while the reference manual states that these should lie between 30 and 
100 degrees. 

The mesh was built in such a way that adjacent elements, like elements 1 
and 2 shown in Figure 5, were able to share common nodes. The top side 
exploded mesh view of the three dimensional model created by ATILA DEPL 
program is shown in Figure 6. 

The types of isoparametric elements used in the mesh generation are 
described in [Ref. 3] and listed in the following Table 2. 



Region 


Element 


Geometry 


Fluid 


PRIS15F 


15-node triangular 
prism 


Radiation Damping 


TRIA06R 


6-node triangular 



Table 2. Isoparametric Elements Used in The Three Dimensional Mesh 
Generation 
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Figure 6. Three Dimensional Spherical Model. Topside View 



15 








16 



IV. RESULTS 



A. EXACT ANALYTICAL RESULTS FOR SCATTERED PRESSURE WAVE 

The linearized homogenous wave equation for the propagation of sound 
in ideal (nonviscous) fluids and the time-independent, lossless Helmholtz 
equation for a pressure wave with time-harmonic dependence p(t,r)=p(r)*e jo)t , at a 
position r=(r,9,cp) and time (t), is given by [Refs. 9, 10] : 



The solution of the Helmoltz Equation (1) for the incident (P jn ) and 
scattered (P sc ) pressure fields in the spherical coordinate system shown in 
Figure 7 is given by Equations 2 and 3: 




( 1 ) 



z 



(r,0.(p) 



x 




y 



Figure 7. The Spherical Coordinates (r,0,cp). 
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( 2 ) 




n = Om = -n 



h[ 1 1) (kr)Y n m (e i <t.) 




n = Om = -n 



•h ( n 2) (ki-)-Y™(e, < |,) 



(3) 



where: 

P nm1 and P nm2 are the amplitudes of the n,m components of the incident 
and scattered waves, respectively, 

h^kr) = the nth order spherical Hankel function of the first kind, 
hn 2) (kr) = the nth order spherical Hankel function of the second kind, 
Y n m (0,<t>) = the spherical harmonic of order n,m, which is related to the 
associated Legendre polynomial by the equation 



Application of boundary conditions on the inner surface of the spherical 
fluid mesh for vanishing of the pressure or its normal derivative, provides 
solutions for the scattered wave for the case of a pressure release or rigid 
surface. 

1. Analytical Expression for the Scattered Pressure from a 
Pressure Release Spherical Surface 

The following equations apply for the n,m component: 




P jn (r 1 0,(p)=Y n m (0,(p)-h n ( 1 )(kr) 



(4) 



P sc (r,0,9)=[-hf 1 1) (ka) / hjf^ka)]- Y n m (0,(p)-h( n 2 )(kr) 



(5) 
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where: 



a = the radius of the spherical shell, 
k = the wavenumber of sound, 

2. Analytical Expression for the Scattered Pressure from a Rigid 
Spherical Surface 



The following equations apply: 



P in (r,e,( P )=Y n m (0,9)-h ( n 1) (kr), (6) 

t t 

P sc ( r ' 0 . ( P)=[-(h { n 1)(ka) ] /(h}, 2 )(ka)]].Y^(e,(p)-hf, 2) (kr), (7) 

where (h^ka))' and (hj, 2) (ka))' are the spatial derivatives of h^(kr) and hj, 2) (kr) 
evaluated at the spherical surface. 

Appendix A contains a C program used to calculate analytically the values 
of the pressure scattered/radiated from a vibrating, spherically symmetric 
surface, for the above boundary conditions [Ref. 11]. 

B. NORMALIZED SCATTERED PRESSURE DEVIATION COMPUTED BY 

ATI LA 

Using the requried ATILA input data file, the pressure scattered by a rigid 
spherical surface due to an incident spherical harmonic wave was computed. 
The results were compared with an exact analytical solution computed using the 
C program described in Appendix A. Accordingly, the limitations of the boundary 
elements were obtained. 

Since the scattered pressure is proportional to the incident pressure and 
has angular dependence, the error in the ATILA results at each finite-element 
node was quantified by normalizing the deviation from the exact value by the 
maximum magnitude of the scattered pressure at the same radius. This 
represents the normalized scattered pressure deviation computed by ATILA. 
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C. RESULTS FOR THE THREE DIMENSIONAL SPHERICAL MODEL 

(MODEL 1) 

The following Figures 8 and 9 depict the results of the calculation for the 
original three-dimensional spherical model (model 1) that had been used in the 
previous investigation [Ref. 2], In this model, the fluid is divided into four layers 
of thickness: 

1. Layer a: 0.25 m 

2. Layer b: 0.25 m 

3. Layer c: 0.5 m 

4. Layer d: 0.5 m 

The radius of the scattering (inner) surface (a) is 0.5m: the radius of the radiation 
boundary is 2.0m. The frequency (f) is 474Hz; the corresponding wavenumber 

of sound in the water is k=— =2m _1 . 

c 

Figure 8, from top to bottom, presents the normalized scattered pressure 
deviation computed by ATILA using multipolar dampers for monopolar and 
dipolar incident pressure fields, versus the dimensionless radius (kr) from the 
center of the structure. Figure 9 presents the normalized scattered pressure 
deviation for the quadrupolar incident pressure field. 

The following Table 3 provides the normalized maximum error in 
percentages in the middle of each fluid layer and the angular location of that 
error for the monopolar, dipolar, and quadrupolar incident fields. 

Nodes located at the poles are points with coordinate spherical angles: 0 
(polar) = 000,180 and 9 (azimuthal) = 000, degrees. Nodes located on the 
equator are points with coordinate spherical angles: 0 (polar) = 090 and 9 
(azimuthal) = 000,090,180,270 degrees. Near the equator, nodes are points with 
coordinate spherical angles: 0 (polar) = 075,105 and 9 (azimuthal) = 
090,180,270 degrees. 
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NORMALIZED SCATTERED PRESSURE DEVIATION 
COMPUTED BY ATILA USING MULTIPOLAR DAMPERS 
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Figure 8. Normalized scattered pressure deviation computed by ATILA versus kr 
for monopolar and dipolar incident fields. 
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NORMALIZED SCATTERED PRESSURE DEVIATION COMPUTED BY 
ATILA USING MULTIPOLAR DAMPERS 
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Figure 9. Normalized scattered pressure deviation computed by ATILA versus kr 
for quadrupolar incident field. 
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Type of 
Field 


Scatterer 

Surface 

ka=1.0 


Center 

of 

Layer a 
kr=1.25 


Center of 
Layer b 
kr=1.75 


Center of 
Layer c 
kr=2.5 


Center of 
Layer d 
kr=3.5 


Maximum 

Normalized 

Error 


kr=1.5 


Monopolar 


2.6 


3.8 


5.7 


7.8 


5.8 


8.6 kr=3.0 


4.6 


m=n=0 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Dipolar 


2.2 


3.0 


5.0 


9.2 


8.5 


11.8 kr=3.0 


3.9 


n=1, m=0 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Dipolar 


1.8 


2.2 


3.6 


5.1 


4.0 


5.7 kr=3.0 


3.0 


n=1, m=+1 


Equator 


Equator 


Equator 


Equator 


Equator 


Equator 


Equator 


Quadrupolar 


0.3 


0.4 


0.8 


2.7 


5.3 


5.8 kr=3.0 


0.6 


n=2, m=0 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Quadrupolar 


0.3 


0.4 


0.5 


0.9 


3.7 


3.6 kr=3.5 


0.5 


n=2, m=±1 


Equator 


Equator 


Equator 


Equator 


Equator 


Equator 


Equator 


Quadrupolar 


0.3 


0.5 


0.9 


1.9 


3.6 


3.6 kr=3.5 


0.5 


n=2, m=±2 


Near 


Near 


Near 


Near 


Near 


Near 


Near 




Equator 


Equator 


Equator 


Equator 


Equator 


Equator 


Equator 



Table 3. Percent normalized scattered pressure deviation and location of the 
maximum error for monopolar, dipolar and quadrupolar incident fields. 

From the above table and Figures 8 and 9, it can be concluded that: 

1. The greatest normalized errors appear for the following type of 

fields: 

a. Monopolar, n=m=0: 8.6% at kr=3.0 

b. Dipolar, n=1, m=0: 11.8% at kr=3.0 

c. Quadrupolar, n=2, m=0: 5.8% at kr=3.0 

2. For the above axisymetrical incident pressure fields the 
corresponding location of the maximum error points is close to the poles and 
exactly on the poles for the scatterer’s surface and on the poles for each layer, 
respectively. 
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3. For the cases, where the maximum error points appear on the 
equator or near the equator (dipolar n=1, m=±1, quadrupolar n=2, m=±1), the 
location of the minimum error points is on the poles. 

4. If the maximum error points are ignored for the monopolar and 
dipolar fields then the maximum normalized error is always less than 4%. 

5. If the maximum error points are ignored for the quadrupolar field 
then the maximum normalized error is always less than 3.5%. 

6. The normalized scattered pressure deviation computed by ATILA 
as it is presented in Figures 8 and 9 appears to increase when moving further 
away from the acoustic center. 

The above results for the monopolar, dipolar and quadrupolar incident 
fields indicate that, for a given value of n, the maximum normalized errors occur 
for the axially symmetric type of fields (i.e., for m = 0). Hence, in investigating 
the influence of the fluid mesh inner radius (ka), the element’s size (kL), and the 
fluid mesh outer radius (kR) on the results, only the axially symmetric incident 
waves were analyzed. 

D. INFLUENCE OF FLUID MESH INNER RADIUS 

In order to evaluate the influence of fluid mesh inner radius on the results, 
models 2 and 3 were used. Recall that model 2 is divided into four equal fluid 
layers of thickness L=0.5m and the scatterer’s radius is 0.5m. Model 3 is divided 
into three equal fluid layers of thickness L=0.5m and scatterer’s radius 1.0m. For 
both cases, the radiation boundary and element’s size satisfy kR=2.5 and 
kL=0.5, with k = 1.0m" 1 and R=2.5m. 

Figure 10 presents the influence of the scatterer’s radius (ka) on the 
normalized scattered pressure deviation, versus kr. The monopolar (a,b), dipolar 
(c,d) and quadrupolar (e,f) fields are shown from top to bottom. On the left side 
the scatterer radius is ka=0.5 and on the right side the scatterer radius is ka=1.0. 
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Figure 10. Normalized scattered pressure deviation computed by ATILA versus 
kr. From top to bottom monopoloar (a,b), dipolar (c,d) and quadrupolar (e,f) 
incident field; on the left: scatterer radius ka = 0.5; on the right: scatterer radius 
ka = 1.0, radiation boundary R and element size L satisfy: kR = 2.5 and kL = 0.5 
for all cases. 
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Table 4 summarizes the results for the maximum error in percent in the 
middle of each layer and the node location for the monopolar, dipolar and 
quadrupolar incident fields when the scatterer radius ka varies and kL and kR do 
not. 



Type of 
Field 


Scatterer 

Surface 


Center of 
Layer a 
kr=0.75 


Center of 
Layer b 
kr=1.25 


Center of 
Layer c 
kr=1.75 


Center of 
Layer d 
kr=2.25 


Maximum 

Normalized 

Error 


kr=1 .5 


Monopolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=m=0 
















ka=1.0 


2.4 


— 


3.5 


5.8 


5.7 


7.0 kr=2.0 


4.6 


ka=0.5 


1.7 


2.9 


4.9 


6.9 


5.9 


7.8 kr=2.0 


5.7 


Dipolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=1, m=0 
















ka=1.0 


1.5 


— 


1.9 


3.3 


3.8 


4.2 kr=2.0 


2.5 


ka=0.5 


0.5 


0.6 


1.8 


3.0 


3.9 


4.0 kr=2.0 


2.2 


Quadrupolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=2, m=0 
















ka=1.0 


0.4 


— 


0.6 


1.4 


2.8 


4.1 kr=2.5 


0.9 


ka=0.5 


0.2 


0.3 


0.5 


1.2 


2.6 


4.1 kr=2.5 


0.7 



Table 4. Percent normalized scattered pressure deviation and location of the 
maximum error for monopolar, dipolar and quadrupolar incident pressure field 
when ka=1.0 and ka=0.5, radiation boundary kR=2.5, element’s size kL=0.5. 

From this table and Figure 10 it is observed that: 

1 . For the monopolar field: 

a. The normalized error at a given value of kR decreases as ka 
is increased from 0.5 to 1.0. 

b. The maximum error decreases as the scatterer radius ka is 
increased from 0.5 (7.8%) to 1.0 (7.0%). 

c. For both cases, the maximum error occurs at kR=2.0. 

d. If the maximum error points (poles) are disregarded, then 
the maximum error is less than or equal to 4%. 
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2 . 



For the dipolar field: 

a. The normalized error at kr = 1.0 increases by a factor of 5, 
from about 0.3% to about 1.5% as the scatterer radius ka is increased from 0.5 
to 1.0. 

b. The normalized error at all other radii is essentially 
unchanged by varying ka. 

c. The maximum error still occurs at kr=2.0 as it occurred for 
the monopolar field. 

d. When the maximum error points (poles) are disregarded 
then the maximum error is less than 3%. 

3. For the quadrupolar field: 

a. The normalized error at a given node shows very little 
dependence on the value of ka. 

b. The maximum error remains the same (4.1%) as the 
scatterer radius ka is increased from 0.5 to 1.0, from the acoustic center 

c. For both cases, if the maximum error points are disregarded 
then the error is less than 2%. 

For the above dipolar and quadrupolar fields, the minimum error appears 
on the equatorial points. Moreover, from Figure 10, very similar deviation curves 
for the monopolar (a,b), dipolar (c,d) and quadrupolar (e,f) fields, for ka=0.5 and 
ka=1.0 are observed. Also, we observe that the maximum normalized error on 
the scatterer’s surface increases as the scatterer’s radius increases. 

E. INFLUENCE OF FLUID MESH ELEMENT SIZE 

In order to evaluate the influence of fluid mesh element size on the 
results, we used models 2 and 4. Recall that model 2 is divided into four equal 
fluid layers of thickness L=0.5m while model 4 is divided into two equal fluid 
layers of thickness L=1.0m. For both cases, the radiation boundary and 
scatterer’s radius satisfy kR=2.5 and ka=0.5, with k = 1.0m _1 , R=2.5m and 
a=0.5m. 
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Figure 11 presents the influence of the element’s size (kL) on the 
normalized scattered pressure deviation, versus kr. The monopolar (a,b), dipolar 
(c,d), and quadrupolar (e,f) fields are shown from top to bottom. On the left side 
the element’s size is kL=0.5 and on the right side the element’s size is kL=1.0. 

Table 5 summarizes the results for the maximum error in percent on the 
scatterer’s surface, in the middle and at the end of each layer and the node 
location for the monopolar, dipolar and quadrupolar fields, when the element’s 
size varies and kR, ka do not. 



Type of 
Field 


Scatterer 

Surface 

ka=0.5 


kr=1.0 


kr=1.5 


kr=2.0 


kr=2.5 


Maximum 

Normalized 

Error 


kr=1.5 


Monopolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=m=0 
















kL=1.0 


1.2 


4.1 


6.3 


6.3 


4.5 


6.3 kr=2.0 


6.3 


kl_=0.5 


1.7 


3.8 


5.7 


7.9 


3.8 


7.9 kr=2.0 


5.7 


Dipolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=1, m=0 
















kl_=1.0 


0.7 


1.0 


2.6 


3.3 


3.9 


3.9 kr=2.5 


2.6 


kL=0.5 


0.6 


0.5 


2.2 


4.0 


3.6 


4.0 kr=2.0 


2.2 


Quadrupolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=2, m=0 
















kL=1.0 


0.2 


0.6 


0.9 


1.7 


4.2 


4.2 kr=2.5 


0.9 


kL=0.5 


0.2 


0.6 


0.7 


1.9 


4.0 


4.0 kr=2.5 


0.7 



Table 5. Percent normalized scattered pressure deviation and location of the 
maximum error, for monopolar, dipolar and quadrupolar incident pressure field, 
when kl_=1.0 and kL=0.5, scatterer radius ka=0.5, radiation boundary kR=2.5. 
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Figure 11. Normalized scattered pressure deviation computed by ATILA versus 
kr. From top to bottom monopolar (a,b), dipolar (c,d) and quadrupolar (e, 
incident field; on the left: element size kL = 0.5: on the right: element size kL 
1.0. Scatterer radius a and radiation boundary R satisfy: ka = 0.5 and kR = 2.5 
for all cases. 
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From this table and Figure 1 1 we observe that: 

1 . For the monopolar field: 

a. On the scatterer’s surface ka=0.5, as the element’s size is 
increased, the maximum error decreases slightly from 1.7% to 1.2%. 

b. The maximum error at each layer becomes greater as kL is 
increased, except at the maximum deviation point, at kr=2.0. 

c. For both cases the maximum error occurs at kr=2.0. 

d. Once the maximum error points (poles) are disregarded, 
then the maximum error is less than or equal to 4%. Also, in that case, the 
corresponding deviation for both cases is similar. 

2. For the dipolar field: 

a. On the scatterer’s surface for both cases the deviation is 
almost the same. 

b. As the element’s size is increased, the error increases 
except at the maximum deviation point at kr=2.0. 

c. The maximum deviation for kL=0.5 occurs at kr=2.0, while 
for kL=1.0 it occurs at kr=2.5. 

d. When the maximum error points (poles) are disregarded, 
then the maximum deviation is very similar and less than 3% for both cases. 

3. For the quadrupolar field: 

a. On the scatterer’s surface, ka=0.5, the maximum deviation 
for both cases is 0.2%. 

b. For both cases, the variation in the normalized deviation 
follows the same trend, and corresponding values are very close to each other. 

c. For both cases, when the maximum error points are 
disregarded, the deviation is less than 2%. 

d. The maximum deviation occurs at kr=2.5 for both cases. 

F. INFLUENCE OF FLUID MESH OUTER RADIUS 

In order to evaluate the influence of the mesh outer radius, we used 
models 2 and 5. Model 2 is the reference model and is divided into four equal 
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fluid layers of thickness L=0.5m, while model 5 is divided into two equal fluid 
layers of thickness L=0.5m. For both cases, the scatterer radius a and elements 
size L, satisfy ka=0.5 and kl_=0.5, for k = 1.0m" 1 and a=0.5m. 

Figure 12, presents the influence of the radiation boundary (kR) on the 
normalized scattered pressure deviation versus kr. The monopolar (a,b), dipolar 
(c,d) and quadrupolar (e,f) fields are shown from top to bottom. On the left side 
the radiation boundary is kR=1.5 and on the right side the radiation boundary is 
kR=2.5. 

Table 6 summarizes the results for the maximum error in percent in the 
middle of each layer and the node location for the monopolar, dipolar and 
quadrupolar fields when the radiation boundary kR varies and kL and ka are the 
same. 

It is observed from this table and Figure 12 that : 

1 . For the monopolar field: 

a. At a given radius, the error increases as the radiation 
boundary kR is increased. For example, in the case where kR=2.5, the error is 
100% greater than the corresponding error when kR=1.5 at the location where 
kr=1.5. 

b. At the outer fluid boundary, the error remains almost the 
same (2%) when the maximum error points (poles) are disregarded. 

c. When kR=1.5, the maximum error occurs at the outer radius 
surface points. On the other hand, when kR=2.5, the maximum error occurs at 
kr=2.0. 

2. For the dipolar field: 

a. The error decreases drastically as we increase the radiation 
boundary kR and becomes almost 100% lower at the boundary for kR=2.5. 

b. The maximum error still occurs at the same radius as it does 
for the monopolar field. 
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Figure 12. Normalized scattered pressure deviation computed by ATILA versus 
kr. From top to bottom monopolar (a,b), dipolar (c,d) and quadrupolar (e,f) 
incident field; on the left: radiation boundary at kR = 1.5; on the right: radiation 
boundary at kR = 2.5; scatterer radius a and element size L satisfy: ka = 0.5 and 
kL = 0.5 for all cases. 
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Type of 
Field 


Scatterer 

Surface 

ka=0.5 


Center of 
Layer a 
kr=0.75 


Center of 
Layer b 
kr=1.25 


Center of 
Layer c 
kr=1.75 


Center of 
Layer d 
kr=2.25 


Maximum 

Normalized 

Error 


kr=1.5 


Monopolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=m=0 
















kR=1.5 


0.8 


1.4 


2.3 


— 


— 


2.3 kr=1.5 


2.3 


kR=2.5 


1.7 


2.9 


4.9 


6.9 


5.9 


7.8 kr=2.0 


5.7 


Dipolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=1, m=0 
















kR=1.5 


0.8 


1.6 


5.3 


— 


— 


8.1 kr=1.5 


8.1 


kR=2.5 


0.5 


0.6 


1.8 


3.0 


3.9 


4.0 kr=2.0 


2.2 


Quadrupolar 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


Poles 


n=2, m=0 
















kR=1.5 


0.2 


0.2 


1.6 


— 


— 


4.8 kr=1 .5 


4.8 


kR=2.5 


0.2 


0.3 


0.5 


1.2 


2.6 


4.1 kr=2.5 


0.7 



Table 6. Percent normalized scattered pressure deviation and location of the 
maximum error for monopolar, dipolar and quadrupolar incident fields, when 
kR=1.5 and kR=2.5, scatterer radius ka=0.5 and element size kL=0.5. 

c. Once again, if we disregard the maximum error points 
(poles), for kR=2.5, the error is always less or equal to 3%. On the other hand, 
for kR=1.5, the error is less than 7%. 

d. The minimum error points for both cases appear on the 

equator. 

3. For the quadrupolar field: 

a. At a given radius, the error decreases as the radiation 
boundary increases. 

b. The maximum error for kR=1.5 and kR=2.5 occurs on the 
outer boundary surface. 

c. For both cases, if the maximum error points (poles) are 
disregarded the error is always less than 2%. 
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d. Finally, the minimum error for both cases appears on the 
equatorial points. 

G. SUMMARY OF RESULTS 

The detailed analyses developed for the five models listed in Chapter III, 
Table 1, indicate that: 

1. For the original three dimensional model (model 1), the maximum 
deviation occurs close to or at the poles for the axially symmetric monopolar 
(8.6%), dipolar (11.8%) and quadrupolar (6%) incident pressure fields. Also, the 
minimum deviation points are at equatorial nodes. Moreover, for the non-axially 
symmetric incident fields, the maximum deviation is less than 5.5% and occurs at 
equatorial points, while the minimum deviation occurs at points located at the 
poles. 

2. Further investigation of the influence of the fluid mesh inner radius 
(ka), the element’s size (kL), and the fluid mesh outer radius (kR), on the results 
for axially symmetric incident pressure fields, revealed that: 

a. For monopolar and dipolar incident fields, when the poles 
are disregarded: 

(1) As the radiation boundary kR increases while ka and 
kL remain constant, the maximum deviation decreases by approximately the 
same amount in percent. This occurs specifically from 6.6% when kR=1.5, to 
4% when kR=2.5, for ka=0.5 and kL=0.5, respectively, which corresponds to 
65%. 

(2) As the scatterer’s radius ka increases, while kL and 
kR remain constant, the maximum deviation decreases by a small amount in 
percent. This occurs, specifically, from 4% when ka=0.5, to 3.7% when ka=1.0, 
for kL=0.5 and kR=2.5, respectively, which corresponds to 8%. 

(3) As the element’s size kL increases, while ka and kR 
remain constant, the maximum deviation remains constant. 

b. For a quadrupolar incident field, the maximum deviation 
remains essentially the same between 1.9% and 1.8%. 
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Table 7, which combines three tables, summarizes our results for the 
maximum normalized error in percent, when the poles are disregarded, and the 
type of field for which this error occurs. Also, from this table it can be observed 
that the maximum deviation occurs for monopolar and dipolar axially symmetric 
incident fields. 



kR=2.5 


kL=0.5 


kl_=1.0 


ka=0.5 


4.0 


4.0 




Monopolar 


Dipolar 


ka=1.0 


3.7 


— 




Monopolar 





ka=0.5 


kL=0.5 


kL=1.0 


kR 


6.6 




1.5 


Dipolar 




kR 


4.0 


4.0 


2.5 


Monopolar 


Dipolar 



kL=0.5 


ka=0.5 


ka=1.0 


kR 


6.6 




1.5 


Dipolar 




kR 


4.0 


3.7 


2.5 


Monopolar 


Monopolar 



Table 7. Percent normalized maximum scattered pressure deviation and type of 
incident pressure field, when poles are disregarded, computed by ATILA. 
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V. CONCLUSIONS AND SUGGESTIONS FOR FURTHER 

INVESTIGATION 



A. CONCLUSIONS 

Five three-dimensional spherical models were developed in order to 
evaluate the radiation boundary elements used in the ATILA finite element code. 
The models simulate a rigid spherical solid structure surrounded by an infinite 
fluid. 

Monopolar, dipolar, and quadrupolar incident spherical pressure wave 
fields were imposed, and the scattered pressure was calculated using ATILA. 
Since the scattered pressure is proportional to the incident pressure and has 
angular dependence, the error in the ATILA results at each finite-element node 
was quantified by normalizing the deviation from the exact value by the 
maximum magnitude of the scattered pressure at the same radius. This 
represents the normalized scattered pressure deviation computed by ATILA. 

The range of values of the dimensionless lengths which characterized the 
problem was: 

The fluid mesh inner radius ka (scatterer radius): 2.0 , 1.0 , 0.5 
The fluid mesh element size kL: 0.5 , 1.0 
The radiation boundary kR: 4.0 , 2.5 , 1.5 

The maximum deviations observed were for the axially symmetric 
monopolar (9%), dipolar (12%), and quadrupolar (6%) incident pressure fields. 
For these cases, it can be concluded that the maximum deviation points are at 
the poles, while the minimum deviation occurs on the equatorial points. 

When the poles are disregarded, for the monopolar and dipolar incident 
pressure fields there is a strong influence of the fluid mesh outer radius (kR) and 
a weak influence of the fluid mesh inner radius (ka) on the results. Specifically, 
as the fluid mesh outer radius increases and/or the fluid mesh inner radius 
increases then the normalized scattered pressure deviation decreases. For the 
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quadrupolar field, the maximum deviation remains essentially constant, 
independent of the above factors (ka, kR). Furthermore, variation of the 
elements size (kL) does not affect the maximum normalized scattered pressure 
deviation. 

Moreover, for the non-axially symmetric incident pressure fields, the 
maximum deviation observed was less than 5.5%, and was found to occur at 
nodes on the equator, while the minimum deviation points were located at the 
poles. 

B. SUGGESTIONS FOR FURTHER INVESTIGATION 

The magnitudes of the errors in the scattered pressures computed by 
ATILA were larger than expected, based upon the mesh sizes employed in the 
models; for all of the models the radial length L of the fluid elements was less 
than 1/6 wavelength. It would have been desirable to examine the calculated 
scattering for a more refined mesh. However, to divide the element scale size by 
a factor of 2 would have resulted in models which exceed the allowed number of 
degrees of freedom. Of course, this is because all the models employed were 
fully three-dimensional, and did not make use of any possible reductions in size 
due to wave function symmetry (in fact, our choice of e'* for the azimuthal wave 
function precluded reduction in that direction). Also, we are interested in the 
computed scattering in all of the spherical harmonic wave components for an 
incident wave of only one component. 

It is suggested that further investigation of the performance of the ATILA 
radiation boundary elements be conducted using two-dimensional models. It 
might also be advantageous to use sin (j> or cos (j> instead of e'™" to represent the 
azimuthal dependence. This will allow a finer mesh to be employed, and a 
broader range of values of ka, kL, and kR to be examined. 
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APPENDIX A. C CODE FOR ANALYTICAL PRESSURE 
SCATTERED CALCULATION 



************* ******* **************************************************** 

* PROGRAM SCHAR.C , by Panagiotis Sinanoglou 2/2/1996 

* Program:Computes the Scattered Spherical Harmonic pressure from a rigid 

* spherical structure at the point r1(x,y,z), for the Wavenumber k, for: 

* natural boundary and pressure release condition 

* Input Variables: 

* x,y,z: Cartesian Coordinates of the Node Point. 

* k: Wavenumber. 

* r2: Radius of the Scatterer's Surface. 

n,m: Orders of Spherical Hankel and Legendre functions. 

* Output Variables: 

* rl: Radius in meters. 

* phi: Azimouthial angle in degrees. 

* theta: Polar angle in degrees. 

* rp5: Real part of computed scattered pressure at rl in pascals 

rp6: Imaginary part of computed scattered pressure at rl in pascals 

* Formula for Incident Wave: Pinc=Pnm(costheta)*e A (imphi)*hn(1)(kr2) 

* Formula for Scattered Wave from Rigid Boundary: 

* Psc=-Pnm(costheta)*e A (imphi)*[hn'(1 )(kr2)/hn'(2)(kr2)]*hn2(kr1 ) 

* Formula for Scattered Wave from Pressure Release Surface: 

* Psc=-Pnm(costheta)*e A (imphi)*[hn(1 )(kr2)/hn(2)(kr2)]*hn2(kr1 ) 

* External Functions: sphbes.c, plgndr.c, bessjy.c, beschb.c, chebev.c, 
complex.c, nrutil.c 

* Header Files for prototype function declaration:nr.h nrutil.h.complexl.h 

* Functions: 

* Normalized Legendre[ plgndr(n,m,x) ], Spherical Hankel for the 

* real and the imaginary parts[ sphbes(n,(k*r2),&xsj,&xsy,&xsjp,&xsyp) ], 

* Complex(structures)[ Cexp(mphi), Cdiv(a,b), Cmul(a,b), RCmul(a.b) ] 

* File nphra33n.c provides node coordinates(x,y,z). 

* File diml.dt stores the corresponding spherical coordinates(r1,theta,phi) 

* File dim13.dt stores the Legendre plgndr(n,m,(z/r1). 

* File dimil. dt stores the real, imaginary part of exp. function e A (imphi). 

* File dim14.dt stores the ratio of the derivatives of the Spherical Hankel for the 
real and 

* imaginary part of hn’(1)(kr2)/hn'(2)(kr2) for natural boundary conditions and the 
ratio of 

*hn(1)(kr2)/hn(2)(kr2) for pressure release, evaluated at the scattering surface, at 
r2=0.5m. 

* File dim15.dt stores the real and the imaginary part of the Spherical 

* Hankel hn2(kr1) at the node range rl. 
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* File dim16.dt stores the real and the imaginary part of the scattered pressure 
Psc, at the * range rl . 

**************** ********************************************* ********** 



#include <stdio.h> 

#include <math.h> 

#include<string.h> 

#include<stdlib.h> 

#include "complexl .h" 

#define NRANSI 
#include "nr.h" 

#include "nrutil.h" 

#define pi 3.141592654 

main() 

{ 

FILE *f1 ; FILE *f3; FILE *f2; FILE *f4; FILE *f5;FILE *f6;FILE *f7; 
FILE *f20; FILE *f21;FILE *f22; 
fcomplex a,b,hn2; 
float r2=0.5,k=2.0; 

float sj,sy,sjp,syp I xsj,xsy,xsjp,xsyp,rp1,rp2,rp3,rp4,rp5,rp6; 
int j.rpm; 

float fac,vaI,t,f11,f12,f13,f14,Pnm,s1,s2; 
unsigned int factorial(unsigned int a); 
double x.y.z.zl ,r,r1 , theta, phi,q,mphi; 
int d,i,N,ch,m,n; 

char g[1 1],h[1 1],c[1 1],e[1 1],f[1 1]; 

f3=fopen("dim1 .dt","w"); 
printf("Enter the integer valuesfor N,n,m\n"); 
scanf("%d %d %d",&N,&n,&m); 
f4=fopen("sphra33b.c","r"); 
f5=fopen("dim1 1.dt”,"w"); 
f6=fopen("dim12.dt","w"); 
f7=fopen("dim13.dt","w"); 
f20=fopen("dim14.dt","w"); 
f21 =fopen("dim 1 5.dt","w"); 
f22=fopen("dim16.dt","w"); 
for(i=0; i<N; ++i) 

{ 

fscanf(f4,"%lf %lf %lf %d 

%lf%s%s%s%s%s\n" 1 &x,&y,&z,&d,&q,&g,&h,&c,&e,&f); 

/* Converts from cartesian(x,y,z) to spherical coordinates(r1, theta, phi).*/ 
if((x>0.0 && y>=0.0)) 

{ r1=sqrt(x*x+y*y+z*z); 
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phi=atan(y/x); 

theta=acos(z/r1);} 

else if(x==0.0 && y==0.0 && z==0.0) 

{r1=0.0; theta=0.0; phi=0.0;} 
else if (x==0.0 && y==0.0 && z!=0.0) 

{ r1=sqrt(x*x+y*y+z*z); 
phi=0.0; 

rl =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1);} 

else if(x==0.0 && y>0.0 &&(z==0.0||z!=0.0)) 

{ phi=pi/2; 

r 1 =sq rt(x*x+y*y+z*z) ; 
theta=acos(z/r1);} 
else if(y==0.0 && x<0.0 ) 

{ phi=pi; 

r 1 =sq rt(x*x+y*y+z*z) ; 
theta=acos(z/r1);} 
else if(y<0.0 && x==0.0) 

{ phi=3*(pi/2); 
rl =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1);} 
else if ((x<0.0 && y>0.0) ) 

(phi=+(pi)+atan(y/x); 
rl =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1); } 
else if (x<0.0 && y<0.0) 

{ phi=pi+atan(y/x); 
rl =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1);} 
else if (x>0.0 && y<0.0) 

{ phi=(2*pi)+atan(y/x); 
rl =sqrt(x*x+y*y+z*z); 
theta=acos(z/r1);} 

theta=theta*( 180/pi); 
phi=phi*360/(2*pi); 

fprintf(f3," %+4.3lf %+6.3lf %+6.3lf\n",r1 1 (theta),(phi)); 

/* The required value of xx=cos(theta) for the Pnm(xx) is: (z/rl) 7 
/* Calculates the Normalized Legendre function "plgndr(n,m,(z/r1) */ 
fprintf(f7,"%4s %4s %10s %24s\n","n'’,"m" ) "x",”plgndr(n,m,x)"); 
rpm=abs(m); 

if( m>=0 ) 

{ f 1 3=factorial(n-m); 
f14=factorial(n+m); 



41 



si =sqrt((2*n+1 )*f1 3/(4*pi*f1 4)); 

Pnm=(plgndr(n,m,(z/r1))*s1); 

fprintf(f7,"%4d %4d %13.6lf %19.6e\n”,n,m,(z/r1),Pnm); 

} 

if(m<0) 

{ fl 3=factorial(n-rpm); 
f 1 4=factorial(n+rpm); 
s2=sqrt((2*n+1)*f1 4/(4*pi*f1 3)); 
f 1 1 =factorial(n-rpm); 
fl 2=factorial(n+rpm); 
r=pow(-1,rpm); 
t=r*(f11/f12); 

Pnm=(plgndr(n,rpm,(z/r1))*t*s2); 
fprintf(f7,"%4d %4d %13.6lf %19.6e\n" 1 n,m,(z/r1),Pnm); 

} 

/* Calculates the e A (imphi) Real/lmagin parts */ 
mphi=(m*phi)*2*(pi/360); 

fprintf(f6,"%lf %lf \n",Cexp(mphi).r,Cexp(mphi).i); 

/* Calculates the (z/r1),phi(degrees) */ 

fprintf(f5,"%+6.3lf %+6.3lf \n M ,(z/r1 ),(phi)); 

/* Calculates the ratio of the Spherical Hankel for the real and the imaginary 
parts */ 

/* of hn'(1)(kr2)/hn'(2)(kr2) on the scattering surface at r2=0.5 for N.B.C , or the 

7 

/* ratio of hn(1)(kr2)/hn(2)(kr2) for pressure release condition 
*/ 



sphbes(n,(k*r2),&xsj,&xsy 1 &xsjp,&xsyp); 
a.r=xsjp; /* a.r=xsj; */ 

a. i=xsyp; /* a.i=xsy; */ 

b. r=xsjp; /* b.r=xsj; */ 

b.i=-xsyp; /* b.i=-xsy; */ 

fprintf(f20,"%30s \n","The real and the imag. of hn’(1)/hn’(2) is:"); 
fprintf(f20,"%f %f \n",Cdiv(a,b).r,Cdiv(a,b).i); 

/* Calculates The real and the imag. of hn2(kr1) */ 
sphbes(n,(k*r1),&xsj,&xsy,&xsjp,&xsyp); 
hn2.r=xsj; 
hn2.i=-xsy; 

fprintf(f21,"%30s \n","The real and the imag. of hn2(r1) is:"); 
fprintf(f21 ,"%f %f \n",hn2.r,hn2.i); 

/* Calculates The real and the imag. parts of the scattered pressure at the */ 

/* radius (rl) :Psc=-Pnm(costheta)*e A (imphi)*[hn'(1)(kr2)/hn'(2)(kr2)]*hn2(kr1)*/ 
rp5=Cmul(RCmul(-Pnm,Cmul(Cdiv(a,b),Cexp(mphi))),hn2).r; 
rp6=Cmul(RCmul(-Pnm,Cmul(Cdiv(a,b),Cexp(mphi))),hn2).i; 
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fprintf(f22,"%f %f\n",rp5,rp6); 

} 

fclose(f4); 

fclose(f3); 

fclose(f5); 

fclose(f6); 

fclose(f7); 

fclose(f20); 

fclose(f21); 

fclose(f22); 

return 0; 

} 

#undef NRANSI 



unsigned int factorial(unsigned int a) 

{ 

if((a==1)||(a==0)) 
return 1; 
else 

{a*=factorial(a-1); 
return a; 

} 

} 

************************************************************************* 



43 














44 



APPENDIX B. FUNCTION INCPRE(X,Y,Z,K) 



FUNCTION INCPRE(X,Y,Z,K) 



* PROGRAM BY ARTHUR LOBO DA COSTA RUIZ 12/23/93 

* MODIFIED BY STEVEN R. BAKER 3/4/96 



* FUNCTION: 

* COMPUTES THE INCIDENT SPHERICAL HARMONIC PRESSURE AT THE 
*POINT (X,Y,Z),FOR THE WAVENUMBER K 

* VARIABLES INPUT: 

* X.Y.Z: CARTESIAN COORDINATES OF THE POINT 

* K: WAVENUMBER 

* VARIABLES OUTPUT: 

* RADIUS: R IN METERS 

* PHI IN DEGREES (HORIZONTAL ANGLE) 

* THETA IN RADIANS (AZIMUTAL ANGLE) 

* INCPRE : PRESSURE AT POINT IN PASCAL 

* 

* FORMULA USED: 

* 

* INCPPRE= H1(N)*P(M,N)*DCMPLX(DCOSD(M*PHI),DSIND(M*PHI)) 

* THIS MAKES AN OUTGOING INCIDENT WAVE FOR E A IWT TIME 

* DEPENDENCE , AS FOR ATILA 



DOUBLE PRECISION K, X,Y,Z,R, PHI/ THETA, KR 
REAL*8 FI ,PMN(-2:2,0:2),LOUT 
INTEGER N,M,NMAX 
COMPLEX*16 INCPRE, H1(0:2),H1OUT 

* N AND M ARE ORDERS FOR HANKEL AND LEGENDRE FUNCTIONS 

N=0 

M=0 

* HI OUT AND LOUT ARE HANKEL AND LEGENDRE OUTPUTS 

* REF TON AND M 

NMAX=2 

* TRANSFORM CARTESIAN COORDINATES (X,Y,Z) 

* INTO SPHERICAL COORDINATES 

* R(RADIUS), PHI(AZIMUTAL ANGLE) AND THETA(POLAR ANGLE) 

R=DSQRT(X*X+Y*Y+Z*Z) 

PHI=DATAN2D(Y,X) 

IF((X.EQ.O).AND.(Y.EQ.O)) PHI=0.0D0 
THETA=DACOS(Z/R) 
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KR=K*R 



NMAX IS THE MAXIMUM NUMBER OF HARMONICS 
CALL HANKEL1(KR,NMAX,H1) 

SUBROUTINE HANKEL1 RETURNS SPHERICAL HANKEL 
FUNCTIONS OF THE FIRST KIND JN + I YN 
CALL LEGNDR(TH ETA, NMAX, PMN) 

LOUT=F1(N,M)*PMN(M,N) 

SUBROUTINE LEGNDR RETURNS ASSOCIATE LEGENDRE FUNTION 
H10UT=H1(N) 

INCPRE=H10UT*L0UT*DCMPLX(DC0SD(M*PHI),DSIND(M*PHI)) 

IF ((R. LE. 0. 501 ).AND.(R.GT. 0.499)) THEN 
PRINT *,"X,Y,Z =",X,Y,Z 
PRINT*, PHI, THETA 
PRINT *,"K,R,KR =",K,R,KR 
PRINT*, "REAL HANKEL1 =",H10UT 
PRINT *, "LEGENDRE =",LOUT 
PRINT *,INCPRE 
PRINT 



* H***************H 



ELSE 

CONTINUE 

ENDIF 

RETURN 

END 



**************************************************************** 



SUBROUTINE HANKEL1(X,NMAX,H1) 
IMPLICIT REAL*8 (A-H.O-Z) 
COMPLEX*16 H 1(0: NMAX) 



C GIVEN THE VARIABLE X, AND THE MAXIMUM ORDER NMAX, 

C THIS ROUTINE GENERATES THE SPHERICAL HANKEL FUNCTION OF 
THE 

C FIRST KIND HI N FOR ALL N FROM 0 TO NMAX (INCLUSIVE) 

C INPUT: 

C X = DOUBLE PREC. VARIABLE (RADIUS) 

C NMAX = INTEGER MAXIMUM ORDER OF BESSEL FUNCTIONS 
DESIRED 
C OUTPUT: 

C H1(N) = ARRAY OF SPHERICAL HANKEL FUNCTIONS H1N(X), WHERE 
C H1N = JN + I YN 

C THIS ROUTINE IS BASED ON THE RECURSION FORMULA 
C FROM ABRAMOWITZ & STEGUN: 10.1.10 & 10.1.15, PP.438-9 
C THE F'S ARE THE COEFFICIENTS OF ORDER N & -(N+1), 

C THE FO'S ARE OLD F'S, FOR RECURSION 
IF ( X .LE. 0.0D0 ) THEN 
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1-11(0) = DCMPLX(1 ,0D0,-1 .0D35) 

DO 2 N = 1, NMAX 
H1(N) = CMPLX(0.0D0,-1 0D35) 

2 CONTINUE 
RETURN 
END IF 
SX = DSIN(X) 

CX = DCOS(X) 

XINV = 1.0D0/X 
MIN = -1.0D0 
FN = XINV 
FMN = 0.0D0 
FNO = FMN 
FMNO = FN 
DO 4 N = 0, NMAX 

H1(N) = CMPLX( FN*SX + M1N*FMN*CX, -FN*CX + M1N*FMN*SX ) 

T1 = (2*N+1)*XINV 
T2 = T1*FN - FNO 
FNO = FN 
FN =T2 

T2 = -T1*FMN - FMNO 
FMNO = FMN 
FMN = T2 
MIN = -M1N 
4 CONTINUE 
RETURN 
END 

Q **************************************************************** 

SUBROUTINE LEGNDR(THETA,NMAX,PMN) 

IMPLICIT REAL*8 (A-H.O-Z) 

REAL*8 F1,PMN(-NMAX:NMAX,0:NMAX) 

C GIVEN THE VARIABLE THETA, AND THE MAXIMUM ORDER NMAX, 

C THIS ROUTINE GENERATES THE ASSOC. LEGENDRE FUNCTIONS PMN 
C OF THE ARGUMENT COS(THETA) (THETA MUST BE BETWEEN 0 & PI) 
C FOR ALL N FROM 0 TO NMAX (INCLUSIVE) 

C AND FOR ALL M FROM -N TO N (SOME OTHERS SET TO ZERO) 

C INPUT: 

C THETA = VARIABLE (POLAR ANGLE), MUST BE BETWEEN 0 & PI 
(INCL.) 

C NMAX = INTEGER MAXIMUM ORDER OF LEGENDRE FUNCTIONS 
C DESIRED 
C OUTPUT: 

C PMN = DOUBLE PREC. ARRAY, CONTAINS ASSOC. LEGENDRE FNS 
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C THIS ROUTINE IS BASED ON THE RECURSION FORMULAE 

C FROM ABRAMOWITZ & STEGUN 
X = DCOS(THETA) 

SINTHT = DSIN(THETA) 

IF ( SINTHT .GT. 0. ) THEN 
SININV = 1.0D0/SINTHT 
ELSE 

SININV = 0.0D0 
END IF 

C SET VALUES FOR N = 0, 1 (NMAX MUST BE AT LEAST 1) 

PMN(0,0) = 1.0D0 
PMN(1,0) = 0.0D0 
PMN(-I.O) = 0.0D0 
PMN(0,1) = X 
PMN(1,1) = -SINTHT 
PMN(-1,1) = SINTHT*0.5D0 

C IN LOOP, TNP1 = 2*N+1, TNP2FC = (2*N+2)!, MIN = (-1)**(N+1) 
TNP1 = 1.0D0 
TNP2FC = 2.0D0 
MIN =-1.0D0 
DO 4 N = 1, NMAX-1 
TNP1 = TNP1 + 2.0D0 
TNP2FC = TNP2FC * TNP1 * (TNP1+1) 

MIN = -MIN 
DO 3 M = -N, N 

PMN(M,N+1) = (TNP1*X*PMN(M,N) - (N+M)*PMN(M,N-1))/(N-M+1) 

3 CONTINUE 
PMN(N+1,N) = 0.0D0 
PMN(-N-I.N) = 0.0D0 

PMN(N+1,N+1) = (X*PMN(N,N+1) - TNP1*PMN(N,N)) * SININV 
PMN(-N-1,N+1) = M1N*PMN(N+1,N+1)/TNP2FC 

4 CONTINUE 
RETURN 

END 

it****************************** 

FUNCTION F1(NN,MM) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 FACT1.FACT2.PI 
PI=4.0D0*DATAN(1 .0D0) 

FACT1=1.0D0 
DO 10 N1.NN+MM 
10 FACT1=FACT1*I 
FACT2=1.0D0 
DO 20 1=1 .NN-MM 
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20 FACT2=FACT2*I 

FI = DSQRT((2*NN+1)*FACT2/4.0D0/PI/FACT1 ) 

RETURN 

END 



49 




50 




LIST OF REFERENCES 



1. Scandrett, C. L. and Canright, D. R., “Acoustic Interactions in Arrays of 
Spherical Elastic Shells”, NPS Technical Report NPS-53-90-009, June 1990. 

2. Ruiz, A. L. C., “Calculation of the Transition Matrix for the Scattering of 
Acoustic Waves from a Thin Elastic Spherical Shell Using the ATII_A F.E.C., 
Master’s Thesis, NPS, Monterey, CA, 1994. 

3. Institut Superieur d’Electronique du Nord (ISEN), ATILA Finite Element 
Code for Piezoelectric and Magnetostrictive Transducer Modeling Version 5.03, 
User’s Manual, September 1993. 

4. Davies, A. J., The Finite Element Method: A First Approach, Clarendon 
Press, Oxford, 1980. 

5. Zienkiewicz, O. C. and Taylor, R. L. The Finite Element Method, McGraw- 
Hill Book Company, Kerkshire, 1989. 

6. Fenner, Roger T., Finite Element Methods for Engineers, The MacMillan 
Press LTD, London, 1975. 

7. Schwartz, H. R., Finite Element Methods, Academic Press, London, 1988. 

8. Zienkiewicz, O. C. and Newton, R. F., “Coupled Vibrations of a Structure 
Submerged in a Compressible Fluid,” Proceedings of the International 
Techniques, Stuttgart, 1969. 

9. Kinsler, L. E., Frey, A. R., Coppens, A. B. and Sanders, J. V., 
Fundamentals of Acoustics, John Wiley & Sons, 1982. 

10. Ziomek, L. J., Fundamentals of Acoustic Field Theory and Space Time 
Signal Processing, CRC Press, Inc., 1995. 

11. Press, W. H., Tenkolsky, S. A., Vetterling, W. T., and Flannery, B. P., 
Numerical Recipes in C, Second Edition, Cambridge University Press, 1992. 



51 









52 



INITIAL DISTRIBUTION LIST 



1. Defense Technical Information Center 2 

8725 John J. Kingman Rd., STE 0944 

Ft. Belvoir, VA 22060-6218 

2. Dudley Knox Library, 2 

Naval Postgraduate School 

41 1 Dyer Rd. 

Monterey, CA 93943-5101 

3. Chairman, Department of Physics 1 

Code PH 

Naval Postgraduate School 
Monterey, CA 93943-5100 

4. Prof. Steven R. Baker 4 

Code PH/Ba 

Naval Postgraduate School 
Monterey, CA 93943 

5. Prof. Clyde L. Scandrett 4 

Code MA/Sd 

Naval Postgraduate School 
Monterey, CA 93943 

6. Prof. Oscar B. Wilson 1 

Code PH/WI 

Naval Postgraduate School 
Monterey, CA 93943 

7. Dr. Regis Bossut 1 

Institut Superieur d’Electronique du Nord 

41 , Blvd. Vauban 
59046 Lille, Cedex, France 

8. Hellenic Navy General Staff 1 

Papagou 

Athens, Greece 15669 

9. Embassy of Greece - Naval Attache 1 

2228 Massachusetts Avenue, NW 

Washington, DC 20008 



53 



10. LT Panagiotis A. Sinanoglou, Hellenic Navy 1 

Embassy of Greece - Naval Attache 

2228 Massachusetts Avenue, NW 
Washington, DC 20008 

11. Dr. Joseph Blue 1 

Head 

Naval Undersea Warfare Center 
Underwater Sound Reference Division 
Orlando, FL 32856 

12. Dr. Roger Richards 1 

Code 2131 (Transducer Branch) 

Naval Undersea Warfare Center 
New London Detachment 
New London, CT 06355-1644 



54 



DUDLEY KNOX LIBRARY 

NAVAL POSTGRADUATE SCHOOL 
MONTEREY CA 93943-5101 



DUDLFY KKinv I inn * mi , 

O ^/Di5 UU323869 2 




